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KINEMATICAL  GEODESY  II 


*  Introduction 

In  an  earlier  report  (Moritz,  1967)  we  h_ve  investigated  theoretical  questions 
related  to  the  use  of  moving  instruments  for  the  measurement  of  elements  of  the 
gravitational  field,  such  as  airborne  gravimeters  or  gradiometers .  The  main  problem 

Jj 

is  the  separation  of  gravitation  and  inertia,  the  extraction  of  the  gravitational  "signal" 
form  the  inertial  "noise".  We  have  seen  that  such  a  separation  is  rigorously  possible 

$ 

A 

i;  for  the  second  and  higher  derivatives  of  the  potential,  but  not  for  the  gravity  vector 

r 

itself.  In  the  latter  case,  in  aerial  gravimetry,  one  is,  therefore,  obliged  to  reduce 
the  effect  of  inertial  noise  as  much  as  possible  by  using  some  statistical  filtering.  The 
fact  that  a  statistical  separation  of  gravitation  and  inertia  can  never  be  perfect,  is 
mainly  responsible  for  the  fact  that  the  accuracy  of  aerial  gravimetry  is  considerably 
inferior  to  the  sensitivity  of  the  instruments  themselves . 

In  Part  A  of  the  present  report  we  shall  describe  the  principles  of  a  rigtrous 
separation  of  gravitation  and  inertia  for  the  gravity  vector  itself.  This  can  be  done 
by  simultaneously  measuring  the  first  and  second  derivatives,  that  is,  by  combining 
a  gravimeter,  or  accelerometer,  with  a  gradiometer. 

In  (Moritz,  1967)  we  have  also  studied  how  measured  second-order  gradients 
can  be  used  in  geodesy.  We  have  seen  that  geodetically  important  quantities,  such  as 
deflections  of  the  vertical  or  geoidal  heights,  can  be  derived  from  these  measurements 
either  by  line  integration  (somewhat  similar  as  in  astrogeodetic  leveling)  or  by  global 
integration,  using  formulas  analogous  to  Stokes'  integral. 

Unfortunately,  such  integral  formulas  have  severe  shortcomings,  which  make 
their  practical  application  hardly  feasible: 

1.  Second-order  gradients  are  much  more  irregular  than  gravity  anomalies, 
so  that  interpolation  is  difficult. 

2.  For  a  Stokes'-type  formula  global  coverage  would  be  necessary. 
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3.  The  available  information  is  very  incompletely  used,  because  only  one 
of  the  five  independent  components  of  the  second-order  gradient  tensor  enters  into 
a  Stokes’-type  integral  formula. 

4.  It  is  impossible  to  combine  second-order  gradients  with  other  data  such 
as  first-order  gradients  or  gravity  anomalies,  in  a  simple  and  well  defined  way,  and 
to  adjust  for  measuring  errors . 

Since  the  first  report  on  Kinematical  Geodesy  was  written,  however,  a  least- 
squares  method  for  estimating  the  terrestrial  gravity  field  (least- squares  collocation) 
was  developed  (Krarup,  1968,  1969;  Moritz,  1970a),  which  avoids  these  shortcomings 
and  permits  the  optimal  use  of  heterogeneous  data  for  the  determination  of  the  earth's 
figure  and  its  gravity  field. 

In  Part  B  of  the  present  report  we  shall  first  discuss  proposed  methods  for 
the  geodetic  use  of  gradients  and  then  apply  least-squares  collocation  to  the  determina¬ 
tion  of  the  gravitational  field  from  measured  gradients .  This  includes  also  the 
determination  of  spherical  harmonics  by  satellite  gradiometiy. 


PART  A 

SEPARATION  OF  GRAVITATION  AND  INERTIA 
1 .  First  and  Second  Order  Gradients 


Let  the  gravitational  potential  of  the  earth  be  denoted  by  V.  Then  tfr  first 

partial  derivatives,  or  first-order  gradients,  V  -bV/bx  etc.  form  a  vector 

i  x 

£ 
i 

| 

j> 

I 

f 

i 

?, 

which  is  the  vector  of  gravitational  force  on  a  unit  mass. 

Adding  to  V  the  potential  of  centrifugal  force  of  the  earth’s  rotation,  we  get 
the  gravity  potential  W.  The  vector 


(1-1) 


wx’ 

v 

wy 

=  g  = 

vy 

w 

V, 

z 

Z 

+  centrifugal  force 


(1-2) 


is  the  gravity  vector.  Since  the  centrifugal  force  is  given  by  a  simple  analytical  expres¬ 
sion  and  can  be  considered  as  known,  the  determination  of  gravitation  (1- 1)  is  equivalent 
to  the  determination  of  gravity  (1”2). 


where 


We  have 

it 

g  cos  a , 

4*. 

1! 

g  cos  /3 , 

(1-3) 

II 

g  cos  y  , 

g  = 

\/w 2  +  w2  +  w 2 

V  x  y  z 

(1-4) 
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is  gravity  and  cos  a,  cos  /3,  cos  y  are  the  direction  cosines  defining  the  direction  of  die 
vertical  or  plumb  line. 

The  measurement  of  the  vector  (1-2)  thus  consists  in  the  determination  of  g 
and  of  the  direction  of  the  plumb  line .  Gravity  g  is  measured  by  gravimetry.  The 
terrestrial  technique  for  determining  the  direction  of  the  vertical  is  to  define  it  by 
means  of  a  spirit  level  and  to  refer  it  to  a  global  rectangular  coordinate  system  by 
means  of  astronomical  observations  (latitude  and  longitude). 

If  the  z-axis  is  made  to  coincide  with  the  normal  to  the  reference  ellipsoid  and 
if  the  x  and  y  axes  are  suitably  oriented,  then,  obviously. 


*»■  ‘  T 


(1-5) 


are  nothing  else  than  the  components  of  the  deflection  oi  the  ve  rtical.  Such  an  orientation 
of  the  coordinate  system  can  always  be  achieved  by  a  coordinate  transformation  so  that 
deflections  of  the  vertical  can  be  computed  from  first-order  gradients. 

In  aerial  measurements,  the  direction  of  the  coordinate  axes  is  maintained  by 
gyroscopic  stabilization,  and  the  three  components  of  the  gravity  vector  can  be  measured 
by  three  accelerometers.  The  accelerometer  output  will  be  affected  by  inertial  disturb¬ 
ances,  which  must  be  removed  as  much  as  possible  by  statistical  filtering. 

The  best-known  technique  falling  under  this  general  principle  is  aerial  gravi¬ 
metry,  where  the  vertical  component  of  the  gravity  vector  is  measured  by  a  gravimeter 
or  a  vertical  accelerometer,  which  is  basically  the  same;  cf.  (Szabo  and  Anthony,  1971). 
For  a  suggestion  for  determining  deflections  of  the  vertical  by  a  similar  principle  cf. 
(Bradley,  1970). 

9  9 

The  second-order  derivatives,  or  second-order  gradients,  Vxx  =d  V/dx  “ 
etc.  form  a  symmetric  matrix  or  tensor 


V 

V 

V  „ 

XX 

xy 

xz 

V 

V 

v 

yx 

yy 

yz 

V 

V 

V 

zx 

zy 

zz 

This  tensor  contains  five  independent  quantities:  because  of  symmetry  we  have 


V  =  v  v  =  V  V=V. 

yx  xy’  zx  xz  *  zy  yz 


and  in  free  space,  Laplace's  equation 


V  +  V  +  V  =  0 
xx  yy  zz 


(1-7) 


(1-8) 


holds,  so  that  the  9  components  of  the  tensor  (1-6)  must  satisfy  4  conditions. 

Instruments  measuring  second-order  gradients  are  called  gradiometers .  A 
stationary  instrument  of  this  kind  is  the  well-known  torsion  balance,  cf .  (Mueller,  1964); 
for  recent  developments  in  airborne  and  satellite-borne  instruments  cf.  (Anthony,  1971), 
(Forward,  1971),  (Savet,  1970),  (Trageser,  1970). 

One  of  the  most  interesting  features  of  gradiometer  measurements  is  that  sec¬ 
ond-order  gradients  measured  by  moving  instruments  are  purely  gravitational,  inertial 
disturbances  having  no  effect  on  them  provided  the  coordinate  axes  are  gyros copically 
stabilized.  For  an  investigation  of  this  problem,  covering  also  the  general -relativistic 
aspects,  see  (Moritz,  1967). 

Anomalous  Gradients .  -  It  is  convenient  to  approximate  the  gravity  potential 
W  by  a  given  simple  function  U,  called  normal  gravity  potential  and  representing  the 
gravity  potential  of  an  equipotential  ellipsoid  (Heiskanen  and  Moritz,  1967,  sec.  2-13). 
The  difference 


T  =  W  -  U 


(1-9) 


is  called  disturbing  potential,  or  anomalous  potential . 
side  the  earth,  since  the  centrifugal  part  in  W  and  U 
in  (1-9).  For  the  same  reason  we  may  also  write 


It  is  a  harmonic  function  out- 
are  equal  and,  therefore,  cancel 


T  =  V  -  V  , 


(1-10) 
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denoting  by  V  .he  earth's  gravitational  potential  and  by  V  the  normal  gravitational 
potential,  which  is  the  normal  gravity  potential  U  minus  the  potential  of  the  centrif¬ 
ugal  force. 

We  may  form  first  and  second  order  gradients  of  the  normal  potential  U 
(or  V,  respectively),  and  subtract  them  from  the  measured  gradients.  In  this  way 
we  obtain  the  anomalous  gradients 


- 

■ 

- 

T 

X 

T 

Lxx 

T 

1xy 

T 

lxz 

T 

and 

T 

T 

T 
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yx 

yy 

yz 

T 

T 

T_, 

T 
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zx 

zy 

zz 

_  _ 
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If  the  position  P(x,y,z)  of  the  measuring  instrument  is  determined  with  res¬ 
pect  to  the  center  of  the  earth,  then  U  and  its  derivatives  can  be  determined  at  the 
observation  point  F  itself,  so  that  the  quantities  (1-11)  can  be  directly  formed.  This 
is,  at  least  approximately,  the  case  in  satellite  gradiometry  where  the  satellite  position 
is  determined  by  tracking,  and  in  aerial  gradiometry  if  position  is  determined  by  inertial 
navigation  connected  to  points  whose  geocentric  position  is  known. 

If,  on  the  other  hand,  the  vertical  position  of  the  aircraft  is  determined  by 
measuring  the  height  above  ground,  then  it  is  more  appropriate  to  consider  the  normal 
gradients  as  referring  to  a  point  Q  which  is  situated  on  the  plumb  line  of  P  and  whose 
normal  potential  U  is  the  same  as  the  actual  potential  W  of  P,  that  is,  UQ  =  IVp. 

See  (Heiskanen  and  Moritz,  1967,  pp.  245-246,  figure  6-3).  Consider,  for  instance, 
the  second  vertical  gradient  W  ,  letting  the  z-axis  coincide  with  the  vertical  through 
P.  Then 


Tzz=Wzz<P>  -  Uzz<p>’ 


whereas  now  we  rather  compute  a  quantity 


=  WZ2<P> 


Uzz<Q> 


(1-12) 


(1-13) 
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The  difference  between  the  two  quantities  is 


Uzz<P>  -  4  tUNp  - 


(1-14) 


where 


Np  =  PQ 


is  the  vertical  separation  of  geop  and  spherop  at  P;  cf .  Fig.  6-3,  loc.  cit. 
For  an  estimate,  assume  a  spherical  normal  earth.  Then 


u=— , 

r 


U_  =  - 


kM 


U. 


rrr 


6kM 

.4  ' 


where  k  is  the  gravitational  constant,  M  is  mass  of  the  earth,  and  r  is  the  radius 
vector.  Near  die  surface  of  the  earth  we  have  approximately 


r  =  R  =  6371  km, 
kM 

|Ur|  =  2  =  G  =  980  gals  , 

|Uzzz|  “  |urrr|  *  ^  =  l.S  x  10‘7  mgal/m2  . 
so  that,  for  Np  =  100  m  , 

IT*  -T  I  =  1.5  x  10' 5  mgal/m  =  0.15E0tv6s. 

I  ZZ  ZZ  I 

Since  100  m  is  a  maximum  value  for  Np,  this  difference  is  well  below  the  expected 
aerial  measuring  accuracy  of  about  1  EttvBs  .  Differences  in  other  second-order  grad¬ 
ients  are  even  much  less  because  of  the  near-spherical  symmetry. 

Thus  we  may  probably  safely  put 

^xx  "  Txx’  ’  ’  ’  ’  ^zz  ”  ^zz  (1-!5) 
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in  most  cases. 


As  for  the  first-order  gradients,  die  fact  just  mentioned  affects  vertical  grad¬ 
ients  and  is  responsible  for  the  distinction  between  gravity  anomaly  and  gravity  disturb¬ 
ance  and  for  the  use  of  the  former  in  aerial  gravimetry.  This,  however,  is  well  known 
and  need  not  be  discussed  here,  cf.  (Heiskanen  and  Moritz,  1967,  pp.  245-246). 

Transformation  of  Gradients. -  Let  an  orthogonal  coordinate  transformation  be¬ 
tween  two  rectangular  coordinate  systems  xyz  and  be  given  by 


T 

where  A  is  an  orthogonal  matrix  and  A  is  its  transpose.  Then  the  first-order  grad¬ 
ients  transform  like 


and  the  second-order  gradients  transform  like 


For  a  derivation  cf.  (Moritz,  1967  ,  sec.  1.5). 

Let  us  now  consider  spherical  polar  coordinates  r(radius  vector),  0  (polar  dis¬ 
tance),  and  X  (longitude).  They  are  related  to  rectangular  coordinates  xyz  by 
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x  =  r  sin  9  cos  \  , 
y  =  r  sin  9  sin  X  , 
z  =  r  cos  9  . 


(1-19) 


Here  the  z-axis  is  the  axis  of  rotation  of  the  earth,  and  the  x-axis  passes  through  the 
Greenwich  meridian;  the  origin  is  at  the  earth's  center  of  mass . 

Together  with  this  global  system  xyz,  let  us  introduce  a  local  rectangular 
coordinate  system  as  follows.  The  origin  is  at  the  point  P  whose  coordinates 
are  x,  y,  z  or  r,  0,1.  The  £  -axis  coincides  with  the  radius  vector,  the  £  -axis 
points  north,  the  rj  -axis  points  east. 

The  first-order  gradients  in  this  system  are  given  by 


T  =  — £L2 -  s  — -  t 

77  r  sin  9dX  r  sin  9  X  ’ 


(1-20) 


The  relation  between  the  derivatives  with  respect  to  x,y,z  and  to  r,  9,X  is 
found  as  follows .  By  the  usual  rules  of  partial  differentiation  we  have 

T  =  T  x  +  T  y  +  T  z 
r  xr  y'r  zr 

and  similar  for  T  and  T  .  The  derivatives  x  etc.  are  obtained  from  (1-19),  for 
9  x  r 

instance 

x  =  sin  9  cos  X  . 
r 
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In  this  way  we  find 


T  =  T  sin  6  cos  X  +  T  sinGsinX  +  T  cos  0 , 
r  x  y  z 

T.  ••  T  r  cos  Geos  X+  T  rcosGsinX"  T  rsinG, 
g  x  y  z 

T.  = -T  rsinG  sinX  +  T  rsinGcosX  . 

X  x  y 

Substitution  into  (1-20)  gives 


with 


-  cos  0  cosX 

-  cos  0  sinX 

sinG 

T 

A  =  -  sin  X 

cos  X 

0 

sin  0  cos  X 

sin  0  sin  X 

cos  0 

(1-21) 


(1-22) 


which  determines  the  transformation  matrix  A . 

Hie  first  derivatives  with  respect  to  r,  0,  X  are  expressed  in  terms  of 
rectangular  gradients  by 


T.  =  r  sin  0  T  . 
X  rj 


(1-23) 
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For  the  second  aerivatives  we  have 


T 

rr 

s 

Tcc  • 

T  . 
r0 

= 

'rT«c‘Tc  • 

T  % 
rX 

= 

r  sin  0  T  „  +  sin  0  T  , 

VC  ri 

T 

00 

= 

i  -  rTj,  , 

T 

ex 

s 

2 

-r  sin0Tk  +  r  cos 0 T  , 

tn  n 

T 

XX 

= 

2  2 

r  sin  0  T  +  r  sin  0  cos  0  Tfc 
W  £ 

(1-24) 


.  2 


These  equations  are  readily  derived  by  differentiating  equations  (1-21)  with  regard 
to  r,  0,  X,  using 


T  =  Tx+Ty+Tz,  etc., 
xr  xx  r  xyJi  xz  r 


and  expressiqg  the  derivatives  with  respect  to  x,  y,  z  by  the  derivatives  with  respect  to 
£ ,  T],  C  by  means  of  (1-17)  and  (1-18)  with  (1-22). 

Conversely  we  have 

Tu  =  ?  Te0  +  7  Tr  • 


1 


T  + 


_  1^.1 
T.  „  = - T 


1  "  r?  2 
W  r  shr  9 


r“  sin  0  ®X  r^sin^0  X 

,  1  T 
r0  +  Te  • 

\  .  +  —  T  +  -V 
XX  r  r  r^ 

1 


1 


Tx  x  +  —  T  +  — cot  0  T  , 

0 


1 


T,  , 


t»C  r  sin0  rX  r*sin0  X 


(1-25) 


T  „  =  T 
CC  rr 
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These  relations  are  found  by  solving  (1-24)  with  respect  to  etc.  and  substituting 

(1-20). 

The  dominant  terms  in  (1-24)  and  (1-25)  are  the  first  terms  on  the  right-hand 
side.  The  following  terms,  representing  effects  of  first-order  gradients,  are  usually 
below  1  Eb'tvo's .  For  instance,  take  the  second  term  on  the  right-hand  side  of  the 
first  equation  of  (1-25).  For  r  =  6371  km  and  =  100  mgals  it  amounts  to 

7-Tr  =  6871  "  °-016  n«al/km 

=  0. 16  EOtvos  . 

We  have  given  these  formulas  for  later  reference.  The  meaning  of  the 
coordinate  systems  introduced  is  as  follows.  The  system  xyz  is  customarily  used 
as  global  coordinate  system  in  geodesy.  The  gradients  referred  to  this  system  can 
be  measured  if  the  inertial  platform  carrying  the  accelerometer  or  gradiometer  is 
made  to  main*ain  a  fixed  orientation  in  inertial  space,  the  direction  of  the  instrument 
axes  coinciding  with  the  directions  of  the  xyz  axes. 

Usually,  however,  the  instrument  axes  are  made  to  slowly  rotate  in  such  a 
way  as  to  always  coincide  with  the  normal  to  the  reference  ellipsoid,  the  tangent 
to  the  ellipsoidal  parallel,  respectively.  For  the  quantities  of  the  anomalous  gravity 
field,  such  as  the  disturbing  potential,  gravity  anomalies,  deflections  of  the  vertical, 
or  anomalous  gradients,  the  spherical  approximation  may  be  used,  which  consists  in 
neglecting  the  flattening  in  ellipsoidal  formulas  so  that  formally  spherical  formulas 
are  obtained.  Loosely  speaking  we  may  say  that  the  reference  ellipsoid  is  formally 
replaced  by  a  sphere  (Heiskanen  and  Moritz,  1967,  sec.  2-14).  Then  spherical 
coordinates  r,  9,  X  may  be  conveniently  used;  and  the  instrument  axes  coincide  with 
the  axes  f,  r>»  C  as  defined  above . 

For  a  summary  of  relevant  aspects  of  inertial  navigation  cf .  (Schultz  and 
Winokur,  1969,  pp.  4892-5). 
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Relations  Between  First  and  Second  Order  Gradients.  -  From  the  second 
derivatives  it  is  possible  to  obtain  the  first  derivatives  by  an  integration  along  the 
flight  path.  For  instance, 

Vx  ■  <Vo  +  /  <Vxxdx  +  Vy+ Vxzdz)-  <l’26> 

Po 

Here  (V  )  refers  to  an  initial  point  along  the  flight  path,  and  V  refers  to  a  current 

X  O  X 

point  P.  If  the  flight  path  is  given  as  a  function  of  the  time  t: 

x  =  x(t) ,  y  =  y(t)  ,  z  =  z(t),  (1-27) 


then 


dx  -  xdt ,  etc. , 


the  dot  denoting  differentiation  with  respect  to  time,  and  we  obtain 

P 

(V  x  +  \  . 

xy 


r 

V  =  (V  )  +  f  (V  X  +  V  y  +  V  Z  )dt, 

x  x  o  J  xx  xy  xz 


V  = 


(V  ) 

y 


+  /  (V  x  +  V  y  +  V  z  )dt  , 

o  J  yx  yyJ  yz 


(1-28) 


r 

V  =  (V  )  +  f  (V  x  +  V  y  +  V  Z  )dt 
z  z  o  J  zx  zyJ  zz 

P  n 


In  order  to  apply  these  formulas,  it  is  necessary  to  know  die  flight  path  as  a 
function  of  time .  In  the  next  section  we  shall  see  how  this  can  be  achieved  using  a 
combined  accelerometer-graaiometer  system. 
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2.  Separation  of  Gravitation  and  Inertia  by 
Using  a  Combined  Accelerometer-Gradiometer  System 

Assume  that  all  first-order  and  all  second -order  gradients  are  simultanously 
measured  in  an  airplane,  die  instrument  axes,  now  denoted  by  Xp  X2»  x^,  being  inert- 
ially  stabilized  in  such  a  way  as  to  maintain  a  fixed  orientation  in  space. 

In  this  case,  the  second-order  gradient  tensor  measured  is  purely  gravitational, 
the  inertial  part  being  zero  (Moritz,  1967,  p.  28);  it  is,  therefore,  given  by  (1-6).  The 
measured  first-order  gradient  vector,  however,  is  not  (1-1)  because  it  is  affected  by 
inertial  disturbances . 

By  equation  (67)  of  (Moritz,  1967)  we  t  **ve 


^7  =  <wik*jk  +  wij)xi  +bi 


(2-1) 


Here  all  subscripts  i,j,k  assume  values  1,2,3;  the  Einstein  summation  convention  is  used. 

* 

The  vector  L  is  the  total  measured  force;  bj  is  the  inertial  disturbance,  the  second 
time  derivative  of  the  position  vector  bj  of  the  origin  of  the  local  frame,  hi  the  present 
case,  the  vector  b^  consists  of  the  three  coordinates  of  the  aircraft  (more  precisely, 
of  the  center  of  mass  of  the  measuring  instrument)  in  a  fixed  coordinate  system.  The 
tensor  describes  the  rotation  of  the  local  frame  with  respect  to  the  fixed  coordinate 
system;  since  we  have  assumed  inertial  stabilization,  w^j  is  zero. 

Thus  (2- 1)  reduces  to 

aV 

SFi  =  f  *  +  b4  .  (2-2) 


Let  us  put  for  the  gravitational  gradients 


av 

SXi  =  v4 , 


axj axj 


=  vij 


(2-3) 
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which  is  a  second-order  linear  differential  equation,  or  rather  a  system  of  linear 
differential  equations,  for  the  velocity  u^ 

To  get  a  clear  picture  of  the  basic  equation  (2-10),  we  shall  write  it  explicitly: 


“x  -  <vxx“x  +  Vy  +  W>  +  Fx  =  °- 

“y  •  <VyxUx  +  vyy“y  +  Vz  >  +  Fy  =  °> 

“x  '  <Vzx“x  +  Vzy“y  +  Vzzuz  >  +  K  *  "• 


(2-U) 


where  uv,  u  ,  u_  are  the  velocity  components  and  F.  F  .  F  are  the  components  of 
x  y  z  x  y  z 

the  measured  force. 

The  quantities  Fj  and  V..  being  given  by  measurement,  eq.  (2- 10)  may  be  solved 
by  die  usual  numerical  methods,  for  instance  by  a  Runge-Kutta  procedure,  to  get  u. . 

This  integration  may  be  performed  in  real  time. 

In  this  way  we  have  indeed  effected  a  separation  between  gravitation  and 

inertia: 

Ihe  gravitational  first-order  gradient  V^,  free  from  inertial  noise,  is  obtained 
from  (2-6)  or,  alternatively,  from  (2-8). 

The  inertial  acceleration,  free  from  gravitational  disturbances,  is  obtained 


as 


bi  =  ui  • 


(2-12) 


It  may  be  twice  integrated  with  respect  to  time  to  give  the  position  vector  b.  in  a  global 


Cartesian  coordinate  system.  Alternatively  we  have  simply 

t 

b,  =  (b.)  +  f  u.  dt  . 

i  io  J  1 


(2-13) 
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Thus  our  combined  accelerometer-gradiometer  system  acts  at  the  same  time 
as  a  purely  gravitational  gravimeter  and  a  true  inertial  navigation  system  that  is  not 
affected  by  the  gravitational  field. 

Additional  Remarks.-  In  order  to  make  the  basic  concepts  transparent,  we 
have  introduced  two  simplifications  which  can,  however,  easily  be  taken  into  account. 

First,  with  current  gravimeter  (or  accelerometer)  and  gradiometer  systems, 
the  force  Fj  or  the  quantities  Vjj  are  not  the  direct  output  when  they  are  functions  of 
time.  In  fact,  the  instruments  may  be  considered  as  linear  oscillating  systems,  for 
which  the  output  is  obtained  as  a  linear  operation  on  the  input  (Fj  or  Vij): 

0  =  L<o  ,  (2-14) 

where  $  is  the  output  and  < p  is  the  force  Fj  or  gradient  Vjj  to  be  measured.  If  the  linear 
operator  can  be  inverted,  then  <p  is  obtained  as 

<p  =  L_10  .  (2-15) 

The  form  of  the  operators  L  and  L  *  is  characteristic  of  the  measuring  system  under 
consideration. 

Secondly,  we  have  already  mentioned  in  sec.  1  that  the  directions  of  the 
instrument  axes,  instead  of  being  kept  fixed,  may  be  rotated  in  a  prescribed  way  such 
that,  for  instance,  one  axis  always  coincides  with  the  normal  to  the  reference  ellipsoid. 
Then  the  rotation  tensor  Wjj ,  instead  of  being  zero,  will  be  given,  so  that  its  effect 
can  be  fully  taken  into  account  using  formulas  such  as  (2- 1).  For  a  special  case  cf . 
(Hansen,  1971). 

Finally  we  remark  briefly  on  the  relation  of  the  present  separation  of 
gravitation  and  inertia  to  the  principle  of  equivalence  of  these  forces .  As  we  have  seen 
in  (Moritz,  1967,  p.  56),  it  is  in  fact  impossible  to  separate  gravitation  and  inertia  as 
long  as  the  force  acting  at  a  point  only  is  considered.  As  soon  as  we  have  a  region  in 
space,  even  an  arbitrarily  small  one,  however,  such  a  separation  becomes  feasible. 
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In  the  present  case,  we  are  given  the  force  along  a  line.  This  fact  by  itself 
is  not  yet  sufficient  for  a  sepai  'tion,  which  becomes  only  possible  through  the  additional 
measurement  of  the  second-order  gradients. 


Here  we  have  restricted  ourselves  to  an  approach  through  classical  mechanics 
because  analyses  such  as  that  in  (Moritz,  1967)  show  that,  to  an  extremely  high  accuracy 
it  gives  the  same  result  as  the  more  rigorous  but  also  more  complicated  approach 
through  the  General  Theory  of  Relativity. 


The  Observational  Data.  -  The  measuring  system  under  consideration  gives  as 
output  the  second-order  gradient  tensor 
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V 

V 

V 

XX 

xy 

xz 

V 

V 

V 

yx 

yy 

yz 

V 

V 

V 

zx 

zy 

zz 

the  first-order  gradient  vector 


V 

x 


LVz 


(2-16) 


(2-17) 


and  even  the  potential  V:  by  integrating 


dV 


V  dx  +  V  dy  +  V  dz 
x  y  z 


we  get 


V  =  Vn  + 


l 

/ 


(V  u  + 

X  X 


V  u 

y  y 


V  u  )dz 
z  z 


(2-18) 


since  the  velocity  components  dx/dt  =  u*  etc.  are  known. 


This  presupposes,  of  course,  that  initial  values  (Vx)0  »  (Vy)0  ,  (VZ)Q  and 
V0  at  some  initial  point  PQ  are  given. 

Likewise  we  obtain  the  position  of  the  aircraft, 


x  =  b.  ,  y  =  b0  ,  z  =  b, 


(2-19) 


as  a  function  of  the  time  t,  again  presupposing  suitable  initial  values  that  were  required 
for  the  integration. 

As  we  have  seen  in  sec.  1,  it  is  convenient  to  subtract  from  the  quantities 
Vij’  Vj  and  V  their  normal  values,  corresponding  to  a  normal  gravity  potential  U, 
to  obtain  anomalous  gradients  and  the  anomalous  potential: 


(2-20) 


T 

T 

T 

XX 

xy 

xz 

T 

T 

T 

yx 

yy 

yz 

T 

T 

T 

zx 

zy 

zz 

We  remark  that,  since  position  is  determined  by  inertial  navigation,  the  normal 
values  will  refer  to  the  same  point  P  as  the  measured  values,  so  that  (Tx»  T^,  Tz) 
represents  the  gravity  disturbance  vector:  see  sec.  1  and  (Heiskanen  and  Moritz,  1967, 
pp.  227  and  245-6). 

Which  of  the  quantities  (2-20)  are  used,  will  depend  on  the  geodetic  computation 

method  chosen.  For  instance,  we  might  use  fT  ,  T  ,  T  )  to  compute  deflections  of 

^  x  y  z 

the  vertical,  to  be  used  for  astronomical  leveling. 

This  would,  however,  be  uneconomical  because  the  available  information  is 
not  fully  used.  In  fact  there  are  5  independent  quantities  (2-20)  since  the  second-order 
gradient  tensor  contains  5  independent  quantities  (see  sec.  1)  and  the  gradient  vector 
and  T  are  obtained  by  integration  of  this  tensor. 
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Thus  we  might  use  five  independent  components  of  the  gradient  tensor,  e.g. 


(2-21) 


or  the  diree  components  of  the  gradient  vector  plus  two  independent  components  of  the 
gradient  tensor,  e.g. 


T 

x 


9 


(2-22) 


or  the  potential  plus  four  other  independent  quantities,  e.g. 


T, 


T  . 
xy 


(2-23) 


It  would,  however,  be  uneconomical  to  use  only 

only  T. 


T  ,  T  or  ,  a  fortiori, 
y  z  ’ 


This  fact  imposes  a  strong  requirement  on  the  geodetic  computation  method 
using  these  data:  in  order  to  take  into  account  all  available  information,  it  should  be 
able  to  use  simultaneously  five  independent  quantities,  ard  it  should  use  them  in  such 
a  way  that  the  result  is  the  same  regardless  of  which  system  of  five  independent  quantities, 
e.g.  (2-21)  or  (2-22)  or  (2-23) is  taken  as  input. 


In  sec.  4  we  shall  present  a  method  that  satisfies  these  requirements. 
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PART  B 


GEODETIC  USE  OF  MEASURED  GRADIENTS 
3.  Review  of  Proposed  Methods 

First-order  gradients  are  equivalent  to  gravity  anomalies  (or  gravity  disturbances) 
and  deflections  of  the  vertical,  as  we  have  remarked  in  sec .  1,  so  that  their  geodetic 
use  may,  in  general,  be  reduced  to  problems  familiar  in  physical  geodesy. 

Still,  this  is  not  the  optimal  procedure,  especially  if  all  first-order  gradients  are 
measured  simultaneously.  The  main  reason  is  that  the  familiar  methods  of  physical 
geodesy  use  either  the  gravity  anomaly  or  the  deflection  of  the  vertical,  a  combination  of 
the  two  types  of  data  not  being  directly  possible.  The  simultaneous  use  of  all  three  com¬ 
ponents  raises  new  problems . 

Such  novel  features  are  particularly  prominent  in  the  geodetic  use  of  second-order 
gradients,  which  is  also  rather  more  difficult.  We  shall,  therefore,  limit  ourselves  to 
considering  second-order  gradients. 

Various  methods  for  their  geodetic  use  have  been  proposed  and  discussed,  e.g. 
in  (Moritz,  1967).  We  shall  now  try  to  give  a  brief  evaluation  of  proposed  methods . 

a)  Line  Integration.  -  We  may  integrate  second -order  gradients  along  the  flight 
path  to  obtain  first-order  gradients .  The  basic  formula  is  (1-28)  or,  written  in  a  more 
gradient  -.ensor  are  measured  and  if  the  velocity  components  Uj  are  known,  either  by 
external  measurements  of  the  flight  path  or  by  the  method  described  in  sec.  2.  The  first- 
order  gradients  so  obtained  are  converted  to  gravity  anomalies  (or  gravity  disturbances) 
and  deflections  of  the  vertical,  which  are  used  in  the  conventional  way  (it  is  e?.sy  to  de¬ 
rive  an  integral  formula  analogous  to  Stokes’  integral  but  using  gravity  disturbances  in¬ 
stead  of  gravity  anomalies). 

The  advantage  of  this  method  is  the  reduction  of  the  problem  to  problems  familiar 
in  physical  geodesy.  A  disadvantage  is  that  the  available  information  is  not  completely 
used:  The  three  components  of  the  gradient  vector  are  computed  from  the  five  independent 
components  of  the  gradient  tensor,  so  that  two  independent  elements  are  not  used. 


Furthermore,  as  we  have  seen  above,  a  combination  of  the  data  furnished  by  the  gradient 
vector  is  not  easily  possible. 

b)  Torsion-Balance  Type  Computations .  -  The  torsion  balance  invented  by 
Edtvtfs,  historically  the  first  and  still  the  only  instrument  in  actual  geodetic  use,  is 
measuring,  not  all  components  of  the  gradient  tensor,  but  only  the  quantities 


V  and  V  -  V 
xy  yy  xx 


(3-D 


(with,  possibly,  Vxz  and  in  addition),  the  xy-plane  being  horizontal.  Hie  quant¬ 

ities  (3-1)  may  be  used  to  calculate  deflections  of  the  vertical  by  an  integration  method 
whose  mathematical  structure  is  clarified  in  (Moritz,  1967,  sec.  1.2). 


This  method,  classical  and  relatively  widely  applied,  is  appropriate  to  the 
torsion  balance.  For  instruments  that  measure  all  components,  the  available  informa¬ 
tion  is  only  partially  used.  Furthermore,  in  this  method  the  lines  of  integration  do  not 
coincide  with  the  flight  path,  so  that  problems  of  interpolation  and  vertical  reduction  oc¬ 
cur  similar  to  those  to  be  discussed  for  the  next  method. 


c)  Global  Integration.  -  This  was  investigated  in  (Moritz,  1967,  sec.  1.3). 
The  relevant  formula  is  equation  (32)  of  that  report: 


T  =  s//TrrSl<*>da  '  <3'2> 

a 

This  integral  formula  is  completely  analogous  to  the  well-known  Stokes  formula:  it 
expresses  the  anomalous  potential  T  in  terms  of  the  second-order  vertical  (radial) 
gradient  Trr  just  as  Stokes’  formula  expresses  T  in  terms  of  the  gravity  anomaly 
Ag  .  The  function  Sj($)  is  a  known  function,  R  is  a  mean  radius  of  the  earth,  do 
is  the  element  of  solid  angle,  and  the  integration  is  to  be  extended  over  the  full  solid 
angle  a,  that  is,  over  the  whole  earth's  surface. 

This  condition,  that  the  integration  be  extended  over  the  whole  earth,  is  the 
more  stringent  as  the  effect  of  the  remote  zones  on  the  integration  decreases  even  less 


22 


than  in  Stokes'  or  Vening  Meinesz'  formulas. 

Since  the  gradients  are  measured  only  at  discrete  points  or  aloqg  certain  profiles, 
they  must  be  interpolated  in  between.  Unfortunately,  second-order  gradients  fluctuate 
much  more  rapidly  and  are  more  irregular  than  gravity  anomalies  or  deflections  of  the 
vertical,  so  that  interpolation  becomes  more  difficult  and  more  problematic.  The  best 
that  can  be  done  to  reduce  interpolation  errors  is  to  use  least-squares  prediction  which 
will  give  as  accurate  results  as  the  data  permit. 

Another  problem  arises  in  this  context.  In  (3-2),  all  the  quantities  Trr  should 
refer  to  the  same  level  surface.  Since  it  will  hardly  be  feasible  to  perform  all  the  meas¬ 
urements  at  the  same  level,  they  might  be  made  at  different  levels  and  reduced  to  the 
same  level.  But  because  of  the  greater  irregularity  of  second-order  gradients,  this  re¬ 
duction  is  even  more  problematic  and  less  reliable  than  for  gravity. 

The  main  objection,  however,  is  that  most  information  remains  unused:  five 
quantities  are  measured  and  only  one  quantity,  Trr,  is  used. 

d)  Determination  of  Spherical  Harmonics .  -  The  anomalous  potential  T  can  be 
expressed  as  a  series  of  spherical  harmonics: 

T  =  E  E  )  n+1(anmcosmX  +  0nm  sinmXjP^cos  0)  ,  (3-3) 

n  =  2  m=0  '  r  > 

where  r (radius  vector),  9(polar  distance)  and  X(longitude)  are  the  spherical  coordinates 
already  used  in  sec.  1,  a  is  the  semi -major  axis  of  the  earth,  (cos  9)  are 
Legendre's  functions,  and  o^m  and  0  are  coefficients  to  be  determined. 

The  use  of  spherical  harmonics  is  most  appropriate  with  satellite  gradiometry 
because  the  convergence  of  the  series  (3-3)  is  satisfactory  at  satellite  altitudes  but  is 
not  so  at  flight  elevations  and  a  forteriori  at  the  earth's  surface,  so  that  an  excessive 
number  of  coefficients  would  be  required  to  get  a  good  approximation  to  the  fine  structure 
of  the  gravity  field. 

By  differentiation  we  may  find  the  corresponding  series  for  the  derivatives  Tr, 
Te»  T^,  Tr0,  TrX,  T00r  Tev  Tu  and  then  the  series  for  T^,  T^,  etc.  by 


(1-25)  or,  alternatively,  the  series  for  ,  T  ,  etc.  by  (1-25)  and 
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which  follows  from  (1-18),  the  matrix  A  being  given  by  (1-22).  In  this  way  we  may 

express  all  measured  second-order  gradients  as  series  which  contain  the  coefficients 

a  and  B  ,  for  instance, 
nm  nm 


T 

xx 


fl(r’  0>  X;%i*  ^nm)  ’ 
f2(r.0,  X ;  0^  . 


(3-5) 


Thus,  every  measurement  gives  one  linear  equation  for  the  c^m  and  /3nm  in 
the  form  of  an  infinite  series;  note  that  r,  9,  X  refer  to  the  particular  point  at  which  the 
measurement  is  performed  and  are  assumed  to  be  known. 

To  determine  the  infinitely  many  o^m  and  Bn m  *  a  finite  number  of  measure¬ 
ments  and,  therefore,  of  linear  equations  is  certainly  not  sufficient.  The  conventional 
procedure  in  this  case  is  to  truncate  the  series  at  some  n  =  np ,  such  that  the  number 
of  retained  parameters  o^m  and  fom  is  smaller  than  the  number  of  observations  and 
these  parameters  can  be  determined  by  an  adjustment. 

Such  a  truncation  is,  however,  a  highly  arbitrary  procedure.  In  the  present 
case  this  is  even  more  problematic  than  in  the  usual  determination  of  spherical 
harmonics  from  orbital  analysis,  since  the  magnitude  of  the  terms  in  the  series  (3-5) 
decreases  considerably  less  than,  e.g. ,  in  the  series  (o  -3),  namely  by  a  factor  of  order 
n  .  Truncation  thus  introduces  "aliasing  errors  "  and  increases  the  mutual  dependence 
of  the  resulting  values . 
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Also  the  statistical  meaning  of  the  formal  adjustment  procedure  is  questionable. 

The  a  and  6  may  be  almost  as  irregular  as  the  effect  of  the  measuring  errors, 
nm  nm 

and  almost  as  small.  It  would,  therefore,  be  more  satisfactory  to  have  a  method  that 
takes  this  fact  into  account  in  a  statistically  well-founded  manner. 

Satellite  gradiometry  will  probably  be  able  to  give  harmonics  of  higher  degrees 
than  does  orbital  analysis;  it  seems,  therefore,  proper  to  comhine  these  two  techniques, 
possibly  with  other  techniques  such  as  satellite  altimetry  and  satellite-to-satellite 
tracking. 

This  brief  discussion  of  different  methods  shows  some  features  that  arise 
particularly  in  the  geodetic  use  of  second -order  gradients: 

1.  A  large  amount  of  information  is  obtained  simultaneously  at  the  same 
point:  the  five  independent  components  of  the  gradient  tensor. 

2.  There  are  difficulties  in  the  application  of  conventional  horizontal  inter¬ 
polation  and  vertical  reduction  techniques  owing  to  the  irregular  and  fluctuating 
nature  of  higher  gradients.  To  better  overcome  these  difficulties,  the  additional 
information  just  mentioned  should  be  used  in  an  appropriate  way. 

3.  By  its  very  nature,  gradiometer  data  are  better  suited  to  give  fine  details 
than  to  provide  the  large  features.  They  are,  therefore,  best  combined  with  other 
data.  This  directly  calls  for  a  method  that  is  able  to  combine  heterogeneous  data  in 

a  natural  way. 

The  classical  methods  of  physical  geodesy-astrogeodetic,  gravimetric, 
dynamic  satellite  techniques --are  always  based  on  data  of  a  single  type.  Attempts  at 
combining  them  are  more  or  less  ad  hoc.  This  is  also  true  for  methods  using  gradio¬ 
meter  data  as  discussed  above,  since  they  are  modeled  after  those  classical  methods. 

4.  Statistical  estimation  and  adjustment  techniques  have  never  penetrated 
very  profoundly  into  classical  physical  geodesy.  Adjustment  techniques  and  methods 
of  error  theory  have  not  been  incorporated  there  in  an  entirely  satisfactory  and 
natural  way.  The  same  holds  for  the  above-mentioned  methods,  which  is  particularly 
serious  here  because  random  errors  may  be  comparable  in  magnitude  to  the  quantity 
to  be  measured,  and  systematic  effects  have  to  be  carefully  eliminated. 
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4.  Application  of  Least-Squares  Collocation 


The  analysis  of  the  preceding  section  has  shown  that  none  of  the  methods  described 
there  is  fully  satisfactory  for  the  geodetic  application  of  gradiometer  measurements .  We 
have  also  recognized  some  desiderata  which  a  better  method  would  have  to  satisfy: 

1.  It  should  be  able  to  handle  all  occurring  data  --  first  and  second  order  grad¬ 
ients  and  any  other  data  —  and  combine  them  in  a  natural,  objective  and  optimal  way. 

2.  It  should  be  able  to  handle  discrete  or  profile  data  at  different  elevations  directly, 
without  interpolation  or  vertical  reduction. 

3.  Methods  described  in  the  preceding  section  should  be  suitable  limiting  cases 
of  it.  For  instance,  if  we  assume  that  only  Trr  has  been  measured,  but  that  it  is  given 
without  errors  at  very  many  points  of  a  level  surface,  then  the  new  method  should  give  a 
result  for  T  that  tends,  as  a  limit  for  infinitely  dense  coverage,  to  the  result  of  eq. 

(3-2). 

4.  It  should  give  the  same  results  whether  second -order  gradients  or  quantities 
derived  therefrom  are  used;  cf.  end  of  sec.  2. 

5.  It  should,  in  a  natural  way,  incorporate  least-squares  adjustment  and  give 
statistically  meaningful  accuracy  estimates.  It  should  be  able  to  make  optimal  use  even 
of  "noisy"  data. 

Recently  a  new  method  of  least-squares  estimation  of  the  gravitational  field  (least- 
squares  collocation)  has  been  developed  which  satisfies  these  requirements  (Krarup,  1968, 
1969;  Moritz,  1970a,  b).  It  may  be  described  as  follows. 

Let  n  quantities  of  the  anomalous  gravity  field  be  measured;  the  measurements 
will  be  denoted  by  Xj,  X2»  ....  xQ.  They  might  be  anomalous  gravity  gradients  but  also, 
e.g.,  conventional  gravity  anomalies,  astrogeodetic  deflections  of  the  vertical  or  geoidal 
heights  derived  from  satellite  altimetry.  Denote  by  Sp  (the  "signal")  the  quantity  of 
the  anomalous  gravity  field  that  we  wish  to  compute,  for  instance  a  geoid  height  or  a 
component  of  the  deflection  of  the  vertical.  Then  Sp  is  given  by  the  matrix  equation 
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Here  C-.  is  the  covariance  between  the  observations  xs  and  Xj  ,  and  C™  is  the 
ij  1  j  n 

covariance  between  die  signal  Sp  and  the  observation  (i,  j  =  1,2, . . .  ,u).  These  co“ 
variances  are  basic;  they  carry,  so  to  speak,  the  burden  of  the  mathematical  structure 
of  the  gravity  field.  Therefore,  much  will  have  to  be  said  in  the  sequel,  especially  in  the 
following  section. 

The  above  formula  presupposes  that  the  (suitably  defined)  average  values  of  sp 
and  x.  are  all  zero: 


M(sp)  =  0,  M(xa)  =  0, 


(4-2) 


which  means  that  both  sp  and  x^  must  be  quantities  of  the  anomalous  gravity  field 
(the  systematic,  "average"  part  of  the  gravity  field  being  removed  by  subtracting  the 
normal  gravity  field)  and  that,  in  addition,  x^  must  not  be  affected  by  systematic  errors . 

The  measurements  Xj  can,  however,  contain  the  effect  of  random  errors;  eq. 
(4-1)  is  valid  in  this  case  as  well  as  in  the  case  of  errorless  observations. 

The  formula  (4- 1)  is  optimal  in  the  sense  that  it  determines  sp  in  such  a  way 
that  the  value  so  obtained  is  compatible  with  the  given  observations  x-  and  the  mean 
square  error  of  estimation  is  a  minimum.  This  has  the  following  meaning.  The  n 
given  observations  do  not  determine  the  gravitational  field  completely  since  this  field 
depends  on  infinitely  many  parameters  (e.g.,  the  full  infinite  set  of  spherical  harmonics). 
Therefore,  there  are  infinitely  many  possible  gravitational  fields  that  are  compatible 
with  the  given  measurements.  To  each  of  these  possible  solutions  there  corresponds 
a  mean  square  error  of  estimation,  mp,  and  eq.  (4-1)  singles  out  that  solution  for 
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which  tDp  is  a  minimum. 


The  formulas  for  this  mean  square  error  of  estimation,  mp  and  for  the  error 
covariances  of  any  computed  values  Sp  and  s^ ,  denoted  by  c^q  ,  are  as  follows : 
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(4-4) 


ihese  quantities  are  analogous  to  the  mean  square  error  after  aujustment  and  the  (error) 
covariance  of  adjusted  values  in  least-squares  adjustment. 


Note  that  in  adjustment  computations,  "variance"  and  "covariance"  always  mean 
error  variance  and  error  covariance,  whereas  in  the  present  method  we  have  both  field 
covariances  (e.g.,  Cpj)  and  error  covariances  (e.g..  o^).  More  about  this  will  be 
said  in  the  next  section. 


For  the  derivation  of  all  these  formulas  see  (Moritz,  1970a,  sec.  2).  If  the 
Xj  are  specialized  to  be  errorless  gravity  anomalies,  the  well-known  formulas  for 
gravity  prediction  result;  note  the  formal  identity  of  the  present  equations  (4-1),  (4-3), 
and  (4-4)  with  equations  (7-63),  (7-64),  and  (7-65)  of  (Heiskanen  and  Moritz,  1967,  sec. 
7-6). 
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Systematic  Effects.  -  Especially  in  moving-base  measurements,  systematic 
trends  such  as  instrumental  drifts  or  systematic  navigation  errors,  are  likely  to  occur. 
They  can  also  be  easily  incorporated  in  the  present  model  by  a  method  developed  in 
(Moritz,  196$ sec.  10)  for  the  case  of  aerial  gravimetry. 


If  the  measurements  Xj  are  affected  by  systematic  errors,  they  are  split  up 
into  a  purely  random  quantity  £•  (comprising  both  signal  and  random  error)  and  a  sys¬ 
tematic  part,  also  called  trend: 

m 


+  E 
a  =  l 


(4-5) 


where  the  X^  are  m  systematic  parameters  and  (A.  ^  denotes  a  given  matrix. 

Thus  the  functional  dependence  on  is  assumed  to  be  linear;  if  it  is  ori 
non-linear,  it  is  to  be  linearized  in  the  usual  way  by  means  of  Taylor's  theorem. 


result 
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The  parameters  Xft 

are  determined  by  a  least-squares  adjustment  with  the 
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are  vectors  or  matrices ,  respectively . 

Then  the  trend  is  subtracted  from  the  data  x ^  to  get  the  "centered  data" 
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and  these  £  may  now  be  used  in  (4-1),  in  the  place  of  x^  to  get  again  an  optimal  esti¬ 
mate. 


A  derivation  of  (4-6)  by  least-squares  adjustment  by  parameters  may  be  found  in 
(Moritz,  1969,  sec.  10).  A  more  satisfactory  simultaneous  deduction  of  (4-1)  and  (4-6) 
from  a  unified  minimum  principle  has  been  given  in  (Moritz,  1970b). 


The  basic  equations  (4-1),  (4-3),  and  (4-4)  need  only  be  slightly  modified  when 
systematic  errors  are  present.  In  (4-1)  we  must  replace  x.  by  ^  as  given  by  (4-8) 
as  we  have  just  seen.  In  (4-3)  and  (4-4),  die  matrix 
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is  to  be  replaced  by 

C_1Ci  -  A (AT  C" 1  A)~ 1 AT C" 1  ]  ,  (4-9) 

where  ^  is  the  n  x  n  unit  matrix,  and  A  and  C  are  given  by  (4-7). 

A  derivation  of  (4-9)  will  be  found  in  the  Appendix. 

Properties  of  the  Solution.  -  As  we  have  remarked  above,  the  present  solution 
is  characterized  by  the  fact  that  the  mean  square  error  of  estimation  is  a  minimum. 
This  is  reminiscent  of  an  analogous  property  of  least-squares  adjustment.  In  fact, 
the  present  method  is  a  generalization  of  least-squares  adjustment  for  the  case  that 
there  is  not  only  a  random  "noise”  (measuring  errors)  but  also  a  random  "signal1 
(elements  of  the  anomalous  gravity  field).  Cf.  (Krarup,  1969)  and  (Moritz,  1970b). 

To  distinguish  it  from  ordinary  adjustment,  the  least-squares  estimation  of  the  gravity 
field  is  called  least-squares  collocation. 

As  we  have  already  mentioned,  the  quantities  x  j,  X2»  ...»  xQ  entering  in  (4-1) 
can  be  any  elements  of  the  anomalous  gravity  field,  affected  or  not  by  random  errors. 
Thus,  eq.  (4-1)  is  able  to  handle  and  to  combine  any  measurements  of  gravitational 
field  elements,  not  only  first  and  second  order  gradients .  Applied  to  gravimetrically 
observed  deflections  of  the  vertical  £,  tj  it  would,  eg.,  give  an  optimally  combined 
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astrogeodetic  and  gravimetric  geoid;  cf .  (Moritz,  1970a,  sec.  9). 

Eq.  (4-1)  could  also  be  used  with  third-order  gradients.  The  reason  why  third- 
order  gradients  are  not  dealt  with  explicitly  in  this  report,  is  that  they  are  probably  of 
less  geodetic  usefulness.  But  the  considerations  of  sec.  1  and  the  techniques  of  secs. 

4,  6,  and  7  could  be  readily  applied  to  third  and  higher  order  gradients  as  well. 

Also  the  signal  Sp  can  be  any  desired  element  of  the  anomalous  gravity  fieid. 

The  different  quantities  computed  in  this  way  are  consistent  with  each  other  in  the 
sense  that  they  belong  to  one  and  the  same  gravity  field. 

In  fact,  the  second  and  third  'actor  in  (4-1),  depending  only  on  the  observations 
Xj  and  their  covariances,  are  the  same  for  all  elements  sp  to  be  computed.  Thus  the 
individual  nature  of  the  quantity  Sp  is  expressed  solely  by  the  first  factor,  the  vector 
(Cp^),  and  the  quantities  Sp  will  be  consistent  if  and  only  if  the  covariances  Cp^  are 
compatible.  The  compatibility  of  these  covariances  is  assured  by  computing  them  ac¬ 
cording  to  the  law  of  propagation  of  covariances  to  be  discussed  below. 

For  instance,  let  all  Xj  be  errorless  measurements  of  the  second  vertical  grad¬ 
ient  Tjj  at  various  points  of  a  level  surface,  and  use  formula  (4-1)  to  compute  Trr 
at  every  point  of  this  level  surface;  this  is,  then,  a  pure  case  of  least-squares  interpola¬ 
tion  in  the  usual  sense.  From  the  continuous  global  Trr-field  obtained  in  this  way,  com¬ 
pute  T  at  some  point  of  the  same  level  surface  by  (3-2).  Alternatively,  compute  T 
directly  from  the  measured  values  using  again  (4-1).  The  resulting  value  for  T 
will  be  the  same  in  both  cases  because  the  covariances  entering  in  (4-1)  are  chosen 
so  as  to  ensure  this. 

In  this  way  we  understand  why  conventional  methods  described  in  sec.  3  can,  in 
fact,  be  considered  as  limiting  cases  of  least-squares  collocation  for  idealized  data 
distributions . 

As  another  example,  consider  the  "problem  of  Bjerhammar":  gravity  anomalies 
are  given  at  discrete  points  of  the  relluroid;  for  a  definition  of  the  telluroid  cf.  (Heiskanen 
and  Moritz,  1967,  p.  292).  As  a  limiting  case,  for  continuous  coverage  of  the  whole 
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telluroid  by  gravity  anoraalies,  this  problem  reduces  to  the  "problem  of  Molodensky  ", 
the  well-known  boundary-value  problem  of  physical  geodesy  (ibid,  p.  291).  The 
Bjerhammar  problem  may  again  be  solved  by  (4-1)  (Moritz,  1970a,  sec.  5);  if  the 
gravity  coverage  becomes  denser  and  denser,  this  solution  tends  to  a  solution  of 
Molodensky  *s  problem.  As,  under  certain  assumptions,  the  solution  of  Molodensky ’s 
problem  is  unique,  this  limiting  solution  will  coincide  with  the  usual  solution  of 
Molodensky 's  problem  by  integral  formulas. 

At  first  sight  it  may  be  difficult  to  believe  that  the  simple  matrix  formula  (4-1) 
is  equivalent  to  complicated  pro''  lures  such  as  the  solution  of  Molodensky *s  problem. 

The  reason  is  that  all  covariances  are  based  on  the  same  covariance  function 
K(P,Q )  (see  next  section),  and  that  this  covariance  function  may  be  selected  to  have  a 
relatively  simple  analytical  expression.  Hence,  all  necessary  operations  may  be  per¬ 
formed  analytically  instead  of  numerically.  Furthermore,  starting  from  the  covariance 
function  of  the  potential,  the  covariances  of  all  relevant  quantities  such  as  gravity  anomalies, 
deflections  of  the  vertical,  or  higher  gradients  are  derived  by  differentiations .  These 
are  much  simpler  to  perform  than  the  integral  operations  necessary  when  going  in  the 
opposite  direction  as  in  the  classical  procedures  of  physical  geodesy. 

By  taking  for  the  covariance  function  a  function  that  can  be  analytically  continued 
down  to  sea  level,  all  difficulties  of  analytical  downward  continuation  are  automatically 
avoided;  such  difficulties  beset  conventional  reduction  procedures. 

These  considerations  help  to  understand  why  (4-1)  is  at  the  same  time  a  generali¬ 
zation  of  classical  procedures,  so  to  speak  with  built-in  interpolation  and  vertical  reduct¬ 
ion,  and  an  essential  simplification. 

There  remains  to  be  discussed  why  the  present  method  gives  the  same  results 

with  any  of  the  data  sets  (2-21),  (2-22)  or  (2-23)  or  with  similar  sets.  The  underlying 

fact  is  that  least-squares  collocation  shares  with  least-squares  adjustment  the  property 

of  invariance  with  respect  to  linear  transformations  both  of  the  signal  Sp  and  of  the  data 

x. .  Invariance  with  respect  to  linear  operations  on  field  elements  Sp  is  the- reason  why 

the  method  determines  a  consistent  gravity  field,  as  we  have  seen  above;  and  invariance 

with  respect  to  linear  operations  on  the  data  is  the  reason  for  obtaining  the  same 

results  with  the  different  data  sets  mentioned,  since  (1-28)  and  (2-18)  are  linear  integral 
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operations.  Cf.  also  (Moritz,  1970a,  pp.  12-13). 

The  equations  of  least-squares  collocation  are  directly  suited  for  high-speed 
computation.  The  biggest  computational  problem  involved  is  the  inversion  of  the 
matrix  C  for  a  great  number  of  observations .  For  a  given  set  of  data  xj,  however, 
such  an  inversion  is  to  be  performed  only  once  for  all  quantities  to  be  computed  and 
for  all  accuracy  evaluations,  as  formulas  such  as  (4-1)  and  (4-3)  show. 

5.  Covariances 

As  we  have  just  seen,  the  covariances  have  ©  carry  the  whole  burden  of  the 
mathematical  structure  of  the  problems  under  consideration.  They  need,  therefore,  be 
investigated  more  closely.  This  has  been  done  in  (Moritz,  1970a,  sec.  4);  we  shall 
summarize  the  relevant  results  and  apply  them  ©  the  present  problem  of  the  use  of 
gradients . 

To  ensure  that  all  our  computed  quantities  belong  to  one  and  the  same  gravity 
field,  all  covariances  that  enter  into  our  computations  must  be  derived  from  a  single 
covariance  function,  for  which  we  may  take  the  covariance  function  of  the  anomalous 
potential  T, 

K(P,Q)  =  cov(Tp,  Tq)  =  M(TpTQ)  ,  (5-1) 

defined  as  the  average  product  of  the  T-values  at  two  points  P  and  Q,  the  average 
L  'jog  understood  in  a  suitable  way. 

The  covariance  function  (5-1)  and  the  quantities  derived  therefrom  are  field 
covariances :  they  express  the  statistical  behavior  of  the  anomalous  gravity  field  and 
should  be  carefully  distinguished  from  error  covariances,  which  express  the  statistical 
behavior  of  observational  errors;  onlv  the  latter  are  considered  in  adjustment  compu¬ 
tations.  Cf.  the  remarks  concerning  the  covariance  function  of  the  gravity  anomalies 
in  (Heiskanen  and  Moritz,  1967,  pp.  267-8). 
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The  results  of  die  computations  do  not  depend  strongly  on  the  choice  of  the 
basic  covariance  function (5-1)  (as  long  as  it  is  used  consistently  throughout !),  in  the 
same  way  as  the  results  in  adjustment  computations  do  not  depend  strongly  on  the 
weights  chosen.  It  is,  therefore,  possible  to  take  for  K(P,Q)  an  analytically  simple 
function. 

K  is  a  function  of  two  points  P  and  Q  defined  on  and  outside  of  some  sphere  of 
radius  R  (which  we  may  take  to  represent  sea  level)  that  must  be  harmonic  both  as  a 
function  of  P  and  as  a  function  of  Q: 

ApK(P.Q)  =  0  =  AQK(P,Q),  (5-2) 

where  Ap  means  the  Laplace  operator  applied  at  the  point  P.  This  follows  immediately 
from  the  definition  (5-1).  Furthermore,  the  function  K  is  assumed  to  be  rotationally 
symmetric:  on  the  sea-level  sphere  of  radius  R,  it  depends  only  on  the  spherical 
distance  $  of  P  and  Q.  Thus 


K(P,Q)  =  K(rp,  rQ,0)  , 


(5-3) 


it  is  a  function  of  the  radius  vectors,  rp  and  rq,  of  P  and  Q ,  and  of  the  spherical 
distance  $  between  P  and  Q. 

Such  a  function  has  a  spherical-harmonic  expression  of  the  form 


K(P,Q)  =  S 
n=0 


n  +  1 

Pn(cos  i/j)  , 


(5-4) 


where  the  k  are  coefficients . 
n 

For  example,  we  may  take  =  k^  =  k^  =  0  and 


k 


n 


A 

(n-l)(n-2) 


for  n  g  3  . 


(5-5) 
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With  these  coefficients,  the  series  (5-4)  may  be  summed  so  that  a  closed  expression 
is  obtained: 


where 


2  .  3 


K(P,Q)  =  A| 


^  j  V 

— j  [P2(cos  0)(l+fn sin2^)]- 


/r  2\  2 
•  *&)  • 


+  A 


cos  0  £n-rz  + 
M 


r0  \  N 


f  [3(t^)  cos*'  1  ] 


2  \  2  _L 

2 


*  [*'  2  (ft-)  cos*+  (-fV)  ] 
\PQ  /  \TQ/  J 


(5-6) 


M  =  1  -  L  -  l~0  )  cost/) 


(5-7) 


N  =  1  +  L  -lA- 

P  Q 


cos  0 


and  A  and  rQ  are  suitable  constants.  According  to  (Lauxitzen,  1971),  to  whom  this 
function  is  due,  it  fits  excellently  global  gravity  and  satellite  data  ,  with 


rQ  =  0.9945  R, 
A  =  7.84888, 


(5-8) 


R  being  again  the  mean  radius  of  the  earth . 

A  simpler  function  which  might  also  be  useful  in  appropriate  cases  is 

1 


K(P,Q)  =  B 


/rPrQ  \ 


rPrQ 


cos  $  +  1 


(5-9) 
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mums 


given  by  Krarup  (1969),  with  suitable  constants  B  and  r^. 

The  equivalent,  for  the  plane,  of  the  spherical  expression  (5-9),  is  the  function 

1 

K(P,Q)  =  c[(xQ-Xp)2+  (yQ-yp)2+ (Zp+ZQ+b)2]"^  (5-10) 


with  constants  C  and  b.  It  is  readily  verified  that  this  function  is  harmonic  with 
respect  to  P  and  Q. 

Propagation  of  Covariances.-  The  law  of  propagation  of  covariances  states 
how  the  covariances  between  any  two  elements  of  the  anomalous  gravity  field  are 
derived  from  the  basic  covariance  function  (5-1).  Perhaps  the  easiest  way  to  express 
it  is  verbally  as  follows: 

Let  u  and  v  be  two  quantities  derived  from  T  by  linear  operations .  Then  the 
covariance  between  u  and  v, 


ccv  (u,  v  )  . 


(5-H) 


is  obtained  as  follows.  Apply  to  the  covariance  function  K(P,Q),  considered  as  a 
function  of  Q,  the  operation  that  do  ;ermines  the  quantity  v  from  T.  To  the  result, 
considered  as  a  function  of  P,  apply  the  operation  that  determines  the  quantity  u 
from  T.  The  result  is  cov(u,  v). 

An  example  will  clarify  this  rule.  Let 


Then  u  is  determined  by  successive  partial  differentiation  with  respect  to  x  and  y, 
and  v  is  derived  from  T  by  partial  differentiation  with  respect  to  z. 

Then,  by  the  verbal  rule  just  given,  cov(Txy ,  Tz )  is  found  as  follows .  Apply 
to  the  covariance  function  K(P,Q),  considered  as  a  function  of  Q,  the  operation  that 
determines  v  from  T,  that  is,  partial  differentiation  with  respect  to  z,  obtaining 
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52K(P,Q) 

To  this  result,  considered  as  a  function  of  P,  apply  successive  partial  differentiation 
with  respect  to  x  and  y,  obtaining 

a2  /«_\s  53k(p,q> 

axpayp  l  3zq  )  9xp3ypSZQ 


Thus  the  desired  covariance  is  given  by 
cov(T  ,  T  )  = 


xy 


a3K(P,Q) 
3xp  Sy0  SZq 


(5-13) 


Putting  x  =Xj,  y  =x2,  z  =x3  and  letting  i,  j,  k,  1  take  the  values  1,  2,  3, 
we  obviously  have 

-(T-  if )  ■  ^  • 

(iX  t\  =  aK<F’Q> 

W1  *\?  ' 

It  *!l  \  =  a2K(H,Q) 

\  ’  aXl  *Xi)  ' 

(  a2T  T\_  *2k(p,Q) 

Vax^Xj’  /  axi,paxjtP 

(XL  ir\  =  a2K(P,Q) 

v***’  axj  J  toijPwjiQ  ' 


covi 


cov 


cov 


cov 


<?T  5T  \  _ 


53K(PlQ  ) 


axk/  axi>pax^p3xk>Q 


cov 


cov 


/  yT 

vaxi  axj 

/ax  a2T  \  =  ,  53k(p, q  )  _ 

3x.ax,J  dxi,P3xj,Q^xk,Q  ’ 

l  /a2T  a?T  \  =  _ a?.K(P.  Q) _ 

Vai^aSr’  Sxjsstj;  axijPaxj>paxktQax1>Q 


(5-14) 
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Here  we  have  used  to  denote  coordinates  x,y,z,  as  we  did  in  sec.  2.  Otherwise 
throughout  Part  B  of  the  present  report  xj(  i  =  1, 2, . . . ,  n  )  always  denotes  measure¬ 
ments,  so  that  no  confusion  should  arise. 

These  formulas  give  the  covariances  between  T  and  its  first  and  second  partial 
derivatives.  Extensions  to  higher  derivatives  are  obvious. 

Sometimes  linear  combinations  occur.  For  instance,  the  gravity  anomaly  is  a 
linear  combination  of  T  and  dT/dr: 

Ag  =  -  -  _L  T  (5-15) 

^  dr  r 

(Heiskanen  and  Moritz,  1967,  p.  89).  Another  example  is  represented  by  (1-25). 

Thus  let  us,  for  instance,  find 

cov(Ag,  ), 

being  given  by  (1-25): 

T«  '  ?Tee  +  fTr  •  <5'16> 

We  shall  use  the  rule  for  the  propagation  of  covariances  as  given  above.  Apply 
to  the  covariance  function  K(P,Q ),  considered  as  a  function  of  Q,  the  operation  that 
determines  from  T  by  (5-16),  obtaining 

4  i  1  * 

rQ  aeQ  rQ  arQ 

To  this  result,  considered  as  a  function  of  P,  apply  the  operation  that  determines  jg 
from  T  by  (5-15).  Thus  we  obtain 
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2 


i  s3k  _  2  <y*K 

■q  8rp  ^2  rp  r2  9  e§ 


i  a^c 

rQ  5rparQ 


rP»Q 


arp 


=  cov(  Ag,  T^)  . 


(5-17) 


In  this  way  we  are  in  a  position  to  express  all  covariances  that  occur  in  the 
geodetic  use  of  gradients,  in  terms  of  partial  derivatives  of  the  basic  covariance  func¬ 
tion  K(P,Q). 


Finally  we  consider  briefly  how  these  partial  derivatives  are  evaluated.  If 
K  is  given  as  a  function  of  rectangular  coordinates  x,y,z,  then  the  evaluation  is 
straightforward.  A  fully  worked  out  example  will  be  found  in  sec.  7;  for  another  example 
see  (Moritz,  1970a,  sec.  7). 


If  K  is  given  as  a  function  of  three  variables  rp,  rp,  0  as  in  (5-3X  then  the 
differentiations  must  be  performed  as 

dxp  arp  *p  +  arp  axp  a^  axp  • 


Now 


r2 

rP 


=  x2  +  y2  + 

P  7p  P 


9 


x2  +  y2  +  z2  , 

Q  yQ  Q 


COS  0  = 


/ 


<P  XQ  +  yP  y0 


^  +y  p  +zp 


+  zpzq 


2 


V^yq+ZQ 


(5-19) 


so  that,  by  straightforward  differentiation, 


drp  Xp 

^p  ’  ^ 


axp  sin  0 


=  0, 


XQ 

rPrQ 


(5-20) 
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In  this  way,  all  occurring  differentiations  may  be  performed  without  ary  mathe¬ 
matical  difficulties,  although  the  analytical  work  may  be  laborious. 

Observational  Errors.  -  If  the  observations  Xj  are  errorless,  then  all  covar¬ 
iances  Cp.  and  Cjj  entering  into  the  basic  collocation  formulas  (4-1),  (4-3)  and  (4-4) 
should  be  field  covariances  as  we  have  just  considered. 

If  the  observations  Xj  are  affected  by  random  errors,  then  the  covariances 
Cp^  remain  field  covariances,  whereas  the  covariances  C-  are  now  given  by 

V  c‘j  +  V  (5'21) 

where  C„  are  the  field  covariances  corresponding  to  the  observed  elements,  and  Djj 
are  the  error  covariances  of  the  observational  errors .  In  the  terminology  of  adjustment 
computations,  the  matrix  (Djj)  is  the  variance -covariance  matrix  of  the  observations. 

Hie  simple  relation  (5-21)  presupposes  that  the  errors  are  uncorrelated  to  the 
anomalous  gravity  field.  This  will  be  true  if  the  observations  have  not  yet  been  sub¬ 
jected  to  a  preliminary  collocation,  for  instance,  a  least-squares  filtering.  In  the 
latter  case,  the  covariance  matrix  (C^  )  is  to  be  taken  from  this  preliminary  collocation; 
cf.  (Moritz,  1969,  sec.  9).  This  is  in  complete  analogy  to  least-squares  adjustment. 

6.  Determination  of  Spherical  Harmonics 

The  collocation  method  described  in  sec .  4  may  also  be  used  to  determine 
spherical  harmonics  from  gradiometer  measurements . 

Let  the  spherical  harmonic  expansion  of  the  anomalous  potential  T  again  be 
given  in  the  form  (3-?),  which  we  shall  write  in  terms  of  fully  normalized  spherical 
harmonics  (Heiskanen  and  Moritz,  1967,  sec.  1-14): 

®  n  /a\  n  +  1  r  _ 

T(r,  e,X)  =n|2  ^=0  (?)  [«nm  R<m 

Then  the  "signal"  Sp  in  (4-1)  is  any  coefficient 


(9'W  +  <6-*> 

or  0  ^  ;  lee  us  assume 
nm 
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(6-2) 


Then 


Cr  =  cov(anm’xi>* 

(6-3) 

C..  =  COV  (Xj  ,  Xj )  , 

(6-4) 

Xj  being  again  any  measured  second-order  gradient  (or  any  other  measured  field  element). 

The  computation  of  the  covariances  Cy  has  already  been  considered  in  the  pre- 
ceding  section;  it  remains  to  study  the  covariances  (6-3). 

The  spatial  covariance  function  of  T  may  again  be  expressed  in  the  form  (5-4): 

K(P,Q)  =  E  11  +  1  Pn(cos0)  .  (6-5) 

n  =2  ^  \  rprQ  /  n 

Then  we  have  by  (Moritz,  1970a,  p.  45) 

_  kn 

cov<5wn’5nm)  =  C0T<'W  *  TS+1  • 

C0''<“nm-V  >  *  V  *  0 

if  p^n  or  q^m  or  both,  (6-6) 

cov(5nm»?pq)  =  0  always. 

In  the  report  just  quoted,  these  formulas  have  been  derived  for  the  covariance 
function  of  the  gravity  anomaly.  It  is,  however,  obvious  that  they  are  valid  for  the 
covariance  function  of  the  potential  as  well. 

Any  gradient  is  obtained  by  single  or  multiple  differentiation  of  T  (or  by  a 
linear  combination  of  such  derivatives),  to  be  symbolized  by  DT.  Thus  from  (6-1) 
we  get 
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•  2  n+l  _ 

DT  =  E  E  a 

n  =2  m  =o 


•  <6*7> 


Let  D  T  denote  the  gradient  the  measurement  of  which  is  Xj .  Then 


=  cov<5nm’  DT) 

(  _ 
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=  cov  <  a  ,  E 
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1  nm  p  =  2 
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,P+1  I  X 
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n  +  l  . 

-  a  cov  ( 


a  ,  a  )  D  (  -1p?1i  |  , 

nm ’  nm  V  ^n+1  J  • 


since  all  covariances  between  coefficients  are  zero  except  one,  by  (6-6).  Thus  we  have 


n+l  /R  <9,X) 

a  i-  i  nm 


Since 


Pi  2n  +  l  n 

Sn,(W) 


kn  D 


n+l 


(6-8) 


.n+l 


=  r"(n  +  1)  Pnm(9)  cos  mX 


is  a  simple  function  of  r,  8,  X ,  any  differentiations  with  respect  to  r,  %  X  are  easily 
carried  out,  e.g.. 


j>L 

d99X 


„n  +  l 


-  (n  +  l)  (n+2)  r  ^n  +  3^  p^  Cos  mX  , 
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d9 


sin  mX  , 
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and  the  components  in  rectangular  coordinates  follow  from  equations  such  as  (1-25). 

For  the  determination  of  the  coefficient  we  need  cov(/?nm ,  DT),  which 

is  again  given  by  (6-8),  with  H  (8 ,  X )  replaced  by  2>(0,  X). 

After  these  preparations  we  are  ready  for  the  application  of  the  collocation  formulas 
such  as  (4-1),  (4-3)  and  (4-4)  for  the  derivation  of  spherical  harmonics  from  gradiometer 
measurements . 

Random  measuring  errors  are  automatically  taken  into  account  if  the  covariance 
matrix  (6-4)  is  properly  computed,  in  the  way  outlined  at  the  end  of  the  preceding 
section . 

Systematic  effects  can  also  be  incorporated  into  our  computations  as  discussed 
in  sec.  4. 

The  advantages  of  the  collocation  method  over  the  conventional  procedure  described 
in  sec.  3  (Item  C)  are  as  follows . 

1.  Every  harmonic  is  determined  independently,  without  aliasing  errors,  since 
the  infinite  series  (6-1)  is  not  directly  vsed  and,  consequently,  no  truncation  occurs. 
Convergence  problems  do  not  affect  the  present  solution. 

2.  The  statistical  meaning  of  the  new  procedure  is  transparent:  it  is  an  optimal 
procedure  in  the  sense  that  it  gives  the  most  accurate  results  obtainable  on  the  basis  of 
the  given  data.  The  statistical  behavior  of  the  anomalous  gravity  field  is  properly  taken 
into  account. 

It  is  said  to  be  a  disadvantage  of  spherical  harmonics  in  satellite  geodesy  that 
their  orthogonality  properties  cannot  be  used  as  efficiently  as  it  would  be  desirable  . 

The  collocation  method  takes  full  advantage  of  these  orthogonality  properties,  in  the  form 
(6-6),  to  separate  the  individual  coefficients . 

Combination  with  any  other  observations --from  classical  techniques  such  as  direct¬ 
ion,  range  and  range-rate  observations  or  from  new  techniques  such  as  satellite  altimetry 
or  satellite-to-sateUite  ranging— are  straightforward  because  (4-1)  can  be  used  with 
heterogenous  observations  as  well,  systematic  parts  being  eliminated  as  discussed  in  sec.  4. 
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7.  Use  of  Profile  Measurements 


Id  (Moritz,  1969,  sec.  6)  we  have  discussed  at  length  the  use  of  aerial  gravity 
measurements  along  parallel  profiles.  Since  the  least-squares  estimation  formulas  hold 
for  any  type  of  measurements,  also  of  different  kinds,  the  formulas  given  there  are  also 
valid  for  measurements  of  first  and  second  order  gradients  along  parallel  profiles. 

To  keep  our  problem  simple  we  assume,  as  we  did  in  the  case  of  aerial  gravi¬ 
metry,  that  the  profiles  are  parallel  straight  lines;  they  need  not  be  equally  spaced,  and 
they  may  be  at  different  elevations.  Let  t  be  the  distance  counted  along  the  direction  of 
the  profiles,  such  that  the  lines  t  =  const,  are  straight  lines  perpendicular  to  this 
direction;  cf.  Figure  2  in  (Moritz,  1969,  p.  26).  Denote  by  x^(t)  the  measurement  of 
some  field  element  (in  our  case,  of  some  first  or  second  order  gradient),  recorded 
along  a  profile  as  a  function  of  t.  In  (Moritz,  1969,  sec.  6),  subscripts  such  as  i  or 
j  (i,j  =  1, 2,  ...,n)  have  labeled  the  profiles;  now  they  label  the  different  quantities  meas¬ 
ured.  All  n  measurements  x^(t)  might,  in  principle,  be  performed  along  the  same 
profile;  or  they  might  be  performed  along  different  parallel  profiles:  the  formulas  are 
the  same. 

The  computational  formulas  derived  in  the  previous  report  just  mentioned  may 
be  summarized  as  follows. 

Denote  by 


CAj(t )  =  COV(xi,  Xj  ) 


(7-1) 


the  autocovariance  function  of  the  measurements .  More  precisely,  C..(t)  is  the 
covariance  between  the  value  of  x^  for  the  argument  u  + 1  and  the  value  of  Xj  for 
the  argument  u,  u  being  any  real  number.  Similarly, 


Cp.(t)  =  COV(Sp,  Xj) 


(7-2) 


denotes  the  cross -covariance  function  between  signal  and  measurement;  more  precisely, 
Cpj(t)  is  the  covariance  between  sp(u  +  t)  and  Xj(u). 
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The  assumption  that  these  covariance  functions  do  not  depend  on  u  but  only  on 
the  argument  difference  (u  + 1)  -  u  =  t  means  that  our  measurements  Xj(t)  are 
considered  as  "stationary  stochastic  processes";  cf.  (Meissl,  1970). 

Now  we  form  the  Fourier  transforms  of  the  covariances,  the  spectra 


cc 

s  (o)  =  /  c  (t)e"1Ut  dt, 

J  -00  J 

« 

S  (U)=  /  Cp.(t)e'1(Jtdt. 

—  oo 


(7-3) 


Next  we  compute  the  "system  functions" 


Hpj(o) 


n  -  /-n 

=  E  S  (U)S.'  %), 
i=l  11 


(7-4) 


where  1\u)  are  the  elements  of  the  matrix  inverse  to  the  n  x  n  matrix  with  elements 
Sjj(u).  Applying  the  inverse  Fourier  transformation  we  obtain  the  "weighting  functions" 


-  00 


and  the  optimum  estimate  of  the  signal  sp(t)  is  finally  given  by 


sp(t)  - 


n 

S 

1  =  1 


J  hp.(t  -0£)x.(a)da  . 


(7-6) 


For  the  validity  of  this  method  it  is  essential  that  the  covariances  be  appropriately 
computed.  If  the  measurements  can  be  considered  as  errorless,  then  all  covariances 
are  directly  given  by  the  law  of  propagation  of  covariances  as  discussed  in  sec.  5.  If 
the  measurements  are  affected  by  random  errors,  Cpj(t)  is  again  a  pure  field  covariance, 
whereas  now 


consists  of  die  field  covariance  Cjj(t)  and  the  error  covariance  Djy(t). 


The  error  covariances  for  second-order  gradients  and  for  first-order  gradients 
obtained  by  the  method  of  sec.  2  should  be  very  much  smaller  than  in  the  case  of  aerial 
gravimetry,  because  there  I\j(t)  also  includes  the  inertial  noise  which  is  now  absent. 

An  Example.  -  This  method  will  be  illustrated  by  a  simple  example.  We  assume 
two  parallel  profiles  1  and  2,  the  first  at  elevation  zj,  the  second  at  elevation  Z2 . 
Along  profile  1,  the  second -order  gradient  Tvv  is  measured,  along  profile  2,  the 
first-order  gradient  Tz  is  measured;  these  measurements  are  errorless. 

This  example  differs  from  Example  2  in  (Moritz,  1969,  pp.  35-36)  only  by 
different  observational  data;  the  geometrical  configuration  (Fig.  1)  and  the 


Figure  1 

mathematical  structure  are  the  same.  The  solution  is  represented  by  equations  (7-12) 
and  (7-13)  on  p.  36  of  that  report.  Only  die  covariance  functions  are  different  because 
of  the  different  observational  data;  they  will  be  computed  now. 

As  the  basic  covariance  function,  let  us  take  the  function  (5-10),  with  C  =  1: 


K(A,B)  =  ~ 


(7-8) 


with 

°2  =  V  +  %  '  yA)2  +  <ZA+  ZB+  b)2' 


(7-9) 
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Both  the  example  and  the  covariance  function  have  been  chosen  for  simplicity;  they  are 
obviously  not  very  realistic. 

By  differentiation  we  find: 


cov(T,  Txy) 

a  k  3  , 

=  <*B3yB  '  • 

(7- 10a) 

cov(T,  Tz) 

9K  1  .  .  .  v 

B+b>- 

(7- 10b) 

cov(T  ,T 
'  xy  xy 

,  * _ als _  * 

aXASyAaVyB 

3  15  ,  .2  15  .  .2 

■  7*  7Vxa>  •7frB'yA)  + 

105  2  2 

+  ^  VXA>  <V  yA>  ' 


cover  .T)  = 


a3K 


15 


=  tV-A-jJVV11  ■ 


xy-z-  8xA8yA9zB  B  A'"E  A  A  B 


covct2.tz>  =  “  ’  pr +  y  W b)2- 


(7-10c) 

(7-10d) 

(7-10e) 


Substituting 

A  =  Pj  ,  B  =  P  ; 

v  xa=  *•  yB"  yA=  ai : 

2  2  2  2 
-  t  +  at  +  (Zp  +  Zj  +  b) 


we  obtain  from  (7- 10a) 


(7- 11a) 


47 


.substituting 


A  =  P2  ,  B  =  P  ; 

V  XA =  '•  V  yA  =  "a2  ’ 

4  m  ,2+a2  +<ZP+Z2"b)2 


we  find  from  (7-  10b) 


CP2(t)“-^-^P+Z2+b)  : 


substituting 


*8  ‘  XA"  '•  V  yA  =  °  ’ 

Djj  =  t2  +(2Zj+  b)2 


we  have  from  (7- 10c) 


Cn(t)  =  TT 
Du 


15 1 


'11 


substituting 


A  =  P2  ,  B  =  Pt  ; 

XB *  XA  *  *>  yB  "  yA  =  "  a  ’ 


2  2  2  2 
°12  =  t  +  a  +(2X+  z2+  bf 


(7- lib) 


(7-llc) 


we  find  from  (7-10d) 


Cl2(t)  - 


15at 


12 


^1+  Z2  +  b) 


(7- lid) 
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and  substituting 


I 

fc‘ 

| 

f 


\ 


xb " xa  =  c »  yB'yA  =  0; 

D22  -  C^^Zj+b)2 

we  finally  obtain  from  (7-lOe) 

C22<‘>  =  -  rT  +  A  <2x2  +  »)2  .  <7- He) 

u22  *^22 

For  the  notations  cf.  Figure  1;  note  that  a,  aj,  a2,  being  measured  in  the  xy -plane, 
are  shown  in  true  size,  whereas  the  spatial  distances  Dj,  D2>  Di2  &ie  shown  as 
projected  onto  the  xy-plane . 

Now  we  can  form  the  Fourier  transforms  of  these  covariance  functions  to  get  the 
spectra  Spl(u) ,  S^u) ,  S^w),  S12(w),and  S22(u).  Then  we  find  H  ^  and  Hp2(u) 
by  eq.  (7-13)  of  (Moritz,  1969,  p.  36),  which  are  nothing  else  than  our  present  eq.  (7-4) 
specialized  for  the  example  under  consideration,  and  hp^(t)  and  hp2(t)  by  (7-5).  Finally, 
(7-6)  gives  T  as  the  signal  to  be  computed. 
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APPENDIX 


We  shall  derive  the  modification  of  the  formulas  (4*3)  and  (4*4)  for  error  variances 
and  covariances  of  the  result  when  systematic  effects  are  present,  arriving  at  (4*9). 
Equation  (4- 1)  may  be  written 


yp-  M- 


(A-l) 


with 


CpC'1 


(A-2) 


we  are  using  a  matrix  notation  similar  to  the  notation  in  (Moritz,  1969,  p.  11),  writing 
yp  for  the  estimated  value  of  Sp  to  distinguish  it  from  the  true  value  sp. 

If  systematic  effects  are  present,  then  in  (A-l)  the  observation  is  to  be  re¬ 
placed  by  the  centered  observation 

£  =  *  "  A  X  ,  (A -3) 

so  that 

yp  =  hpfe'AX)  .  (A -4) 

with  hp  again  given  by  (A-2). 

The  error  of  estimation  is  then  the  difference  true  minus  estimated  value: 

CP"  SP‘  yP  ’ 

and  by  (A -4), 

*p=  SP  ’  llp(iL  *  A  X)  .  (A-5) 

Let  us  now  introduce  the  true  values  of  the  parameters,  X',  and  the  corresponding 
true  values  of  the  centered  observations,  £ ',  for  which  we  have 

(A-6) 


=  x  -  AX'  , 
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in  analogy  to  (A-3).  Substituting 


x  =X'  +  ax' 

(A-7) 

in  (A-5)  we  find 

*p=  sp-  Hpi'-HpAtX'-X)  . 

(A-8) 

The  estimated  values  of  the  parameters  are  given  by  (4-6), 

which  may  be  abb  rev - 

iated  as 

X  =  Hx, 

(A-9) 

T  —  !  -i  T  —  1 

H  =  (ATC  A)  1  A1  C  . 

(A- 10) 

Thus  by  (A-7), 

X  =  Hx  =  Hi'  +HAX', 

and  by  (A- 10)  , 

HAM 

(A-ll) 

(^denotes  again  the  unit  matrix),  so  that 

X'  -  X  =  -  H  j'  , 

(A- 12) 

which  is  substituted  into  (A -8)  to  give 

fp  =  sp  -  hp  (I  -  A  H)  £ '  . 

(A- 13) 

With  the  abbreviation 

EP  =  Up  (I  -  A  H) 

(A- 14) 

we  may  write  this  as 
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Thus 


(A- 15) 


€p  fq  =  Sp  Sq 


h.p  4 


Pi  SQ 


.  h  t't' T  .  T 

+  !lpl£  Sq, 


and  on  forming  the  mean  value: 

°PQ  =  CPQ  '  £pEq  Ep£q  +  EpChq  • 

Now 

Ep  =  hp  (I  -  AH)  =  CpC  ’ 1  [i  -  A ( AT  C  ' 1  A )‘ 1  AT  C“  X] 


(A- 16) 


is  substituted  into  (A- 16)  to  give,  after  some  straightforward  manipulations, 

c^Q=CIq-CpC‘1[l.-A(ATC‘1Ar1ATC'1]cQ  ,  (A-17) 

which  is  (4-4)  with  C -1  replaced  by  (4-9);  and  setting  Q  =  P  gives  the  corresponding 

2 

result  for  the  error  variance  mp . 
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